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Study of the charge correlation function in one-dimensional Hubbard heterostructures 
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We study inhomogeneous one-dimensional Hubbard systems using the density matrix renormal- 
ization group method. Different heterostructures are investigated whose configuration is modeled 
varying parameters like the on-site Coulomb potential and introducing local confining potentials. 
We investigate their Luttinger liquid properties through the parameter K p , which characterizes the 
decay of the density-density correlation function at large distances. Our main goal is the investiga- 
tion of possible realization of engineered materials and the ability to manipulate physical properties 
by choosing an appropriate spatial and/or chemical modulation. 



I. INTRODUCTION 

A key aspect of materials research is to find param- 
eters to tune the physical characteristics of the system 
like conductivity and other desired properties. In the 
last decades there has been enormous progress in the 
generation of nanoscopic quasi-one-dimensional systems, 
e. g., carbon nanotubes |l|, |2j, semiconducting quantum 
wires @, 0] and organic molecules [3J; as well as an in- 
tense study of their transport properties [g, 0, Q such 
as superconductivity Q and quantum Hall edge states 
[lOl Qj| . While the properties of homogeneous one di- 
mensional systems (even with disorder) are relatively well 
understood, very little is known about the properties of 
strongly interacting inhomogeneous systems. Despite of 
the large effort in the study of heterostructures and quan- 
tum dots Qjl, [La, G3, Ha, [l6| , there are still open ques- 
tions which become relevant in modeling the transport 
through molecules, where the electrons interact strongly 
due to the reduced dimension. In addition, its chem- 
istry induces potential barriers which alter the transport 
properties drastically. Technically it is very important to 
know how to control the transport and equilibrium prop- 
erties. In this paper we present a detailed investigation of 
correlation effects in an inhomogeneous one-dimensional 
system including potential barriers. 

The strong electron correlations, inherent to the low- 
dimensional structure, and the large quantum fluctua- 
tions induce new and interesting quantum phases. The 
relevant degrees of freedom are no longer the single par- 
ticle electronic states but the collective spin and charge 
density waves. The low-energy electronic single parti- 
cle excitations possesses vanishing spectral weight at the 
Fermi surface. The physics of such systems, in the ho- 
mogeneous low-energy regime, is well described by the 
Tomonaga-Luttinger liquid (TLL) model [13, HH intro- 
duced by Haldane |19j. Within this model, it is found 
that all correlation functions exhibit a power-law decay 
with the distance, which is specified only by the parame- 
ter K p , known as the Tomonaga-Luttinger (TL) param- 
eter. 

For inhomogeneous structures the high-energy physics 
is determined by the underlying chemistry which, in the 
atomic scale, introduces Coulomb correlations and local 
potentials. On the other hand, at large length scales, 



the physics has to be described by the TLL model. In 
order to establish a connection between the low-energy 
TLL and the quasi-one-dimensional systems synthesized 
in the laboratory, we investigate the density-density cor- 
relation functions in the asymptotic region (i. e. for well 
separated positions x and x'). Position dependent on-site 
Coulomb interaction U(x) and a local potential V(x) are 
used to model the changes in the local chemistry of the 
heterostructures. This defines regions which, for slowly 
varying potentials, can be separately considered as ho- 
mogeneous. We wish to study how the TL parameter 
changes close to the crossover regions. We expect to find 
a description of it in terms of U(x) and the local density 
n(x). 

The paper is organized as follows: In Sec. [IT] the com- 
position of the investigated heterostructures is described 
and we plot our expectations in terms of the coupling 
parameters. In Sec. lIIII we briefly recall the approximate 
results in the low-energy regime for correlation functions 
in the homogeneous case, and we describe the numerical 
procedure, the DMRG method, used to study the one- 
dimensional heterostructures. The results are presented 
and discussed in Sec. IIVI Finally we state our conclu- 
sions. 



II. HUBBARD HETEROSTRUCTURES 

The Hubbard heterostructures we investigate are 
chains with a length of L sites and on-site Coulomb inter- 
action U, which switches between two different values. In 
our case it can be visualized as a valley around the mid- 
dle of the chain with sharp edges at the sites labeled xj, 
and xr . We will refer to this system as Heterostructure 
I. We expect that the slight discontinuity in the charge 
distribution, caused by this form of interaction, will not 
strongly affect the correlation between the adjacent re- 
gions and will make it possible to find a TLL behavior, 
even in the region after the change in the U interaction. 
On a second heterostructure (called from here on Het- 
erostructure IT), in addition to the Coulomb interaction 
described, two potential walls are introduced through the 
confining potential V(^ U). Because of the sharp dis- 
continuity in the charge distribution, we do not expect to 
find a TLL extending beyond the point of the scattering 



potential, however it might still be possible to approxi- 
mate the TLL in the different subchains, since in each of 
them we expect to find a homogeneous particle distribu- 
tion. Fig. Q] shows the layout of the heterostructures. 
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FIG. 1: General arrangement of a Hubbard heterostructure. 
Fhe measurements for (n(x)n(xo)) were carried out from the 
middle point xq — 120. 



Our starting point is an inhomogeneous form of the 
Hubbard Hamiltonian: 



H 



where 
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, a (ci, CT ) is the creation (annihilation) operator 

with spin a (=f, |) on the site i and n icr = c] a Ci^ G is the 
electron number operator, t = 1 is the nearest neighbor 
hopping matrix, which we choose to set the energy scale. 
The Hamiltonian in Eq. |T]) incorporates the different 
systems we want to study and will allow us to find out 
if such systems resemble a TLL and, in that case, also 
to determine the K p parameter from the density-density 
correlation function. 

The sites xl and xr divide the whole system into 
three homogeneous subchains Ul,Uc and Ur, raising 
two questions: first, how the charge correlation function 
behaves in the whole system and second, whether the 
known results for the homogeneous regime can be recov- 
ered within the subchains. 



A. Homogeneous regime 



In a homogeneous TLL, K p determines the long- 
distance decay behavior of all the correlation func- 
tions. In the absence of external magnetic field or spin 
anisotropic interactions, the charge correlation function 
is given by 
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Even though the constant coefficients Ai,A 2 , and B\ 
depend on the model, the algebraic decay is characterized 
only by K p . Of special physical interest are the charge 
density waves with wave vectors 2/cp and 4fcp. While the 
2kp mode dominates over the 4/cp for K p > ~, for suf- 
ficiently large values of the on-site Coulomb interaction 
U, the 4fc F charge mode dominates over the 2fc F mode. 

As a test for our numerics, we considered the case of 
a homogeneous chain for which we confirmed the results 
obtained from the Bethe Ansatz [20, [22[ for the correla- 
tion functions. In Fig. [2] we show our results for several 
values of U obtained with a homogeneous chain of length 
L = 240 sites. We will use this form of the density corre- 
lation function to analyze the low-energy behavior of the 
Hubbard heterostructures. 




III. APPROXIMATE DESCRIPTION OF 

ONE-DIMENSIONAL SYSTEMS IN THE 

LOW-LYING ENERGY SECTOR 



FIG. 2: Tomonaga-Luttinger parameter values for the Hub- 
bard lattice compared to our numerical evaluation of K p 
(dots). U = 1.0, 2.0, 4.0, 8.0, 16.0 from top to bottom. 



The low-lying energy, long-distance physics of one- 
dimensional fermionic systems is described by bosonic 
collective excitations. This bosonization technique yields 
an exact solution for the TL model, whose complete de- 
scription depends solely on the charge and spin velocities 
and the TL parameter K p . In the next subsections we 
first briefly recall the known results [2fJ, [21( for the den- 
sity correlation function in the case of homogeneous sys- 
tems and explain in detail the numerical method used to 
measure the correlation functions in the inhomogeneous 
systems. 



B. Inhomogeneous regime 
Numerical study 

The measurement of observables, which include ground 
state energies and correlation functions, is carried out us- 
ing the density matrix renormalization group (DMRG) 
[23l . |241 . [25( 1 . a method whose roots go back to the nu- 
merical renormalization group formulated by Wilson [26| . 
The DMRG is an efficient numerical method developed 



to overcome the intrinsic difficulties of low-dimensional 
strongly interacting systems. 

The DMRG provides two algorithms to handle an oth- 
erwise exponentially-increasing Hilbert space of a many- 
body system. Both implementations, finite-size and 
infinite-size DMRG base, as in Wilson's renormalization 
group, on a blocking treatment of a lattice system in real- 
space, whose basis of the corresponding Hilbert space is 
decimated under a certain criterion. In the renormal- 
ization group procedure, the decimation of the system's 
basis is done by selecting m states with the lowest energy 
eigenvalues to obtain the ground state of a system. This 
proved to be a reliable method to solve systems, such as 
the Kondo problem, for which the coupling between suc- 
cessive sites decreases exponentially. Thus, it was plausi- 
ble to ignore the connections between neighboring blocks, 
setting to the wave function at the sites outside of the 
block of interest. This lead to inaccuracies when study- 
ing systems such as the Hubbard model, where there is 
no intrinsic separation of the energy scales. To solve this, 
White proposed other criteria to handle both the bound- 
ary conditions when adding a new site to the system as 
well as the selection of states to best represent it. 

The DMRG method considers the system to be con- 
nected to a bath, which is a second block, forming in 
total a superblock. The interactions between the system 
and the bath set the boundary conditions at the edge 
sites of the system as if it would be part of a larger sys- 
tem. In this way, the procedure becomes more accurate 
as the system gets larger. The wave function in the su- 
perblock has the form |\I>) = Ylij^ij\i) ® \J)> where i 
are the states on the system and j are those on the bath. 
From this, the reduced density matrix of the system is 
Pa 1 = ^2j ^ij^i'j- The crucial point is that the density 
matrix contains all the information needed to calculate 
any property of the system and so, the state of the sys- 
tem can be optimally represented by keeping the m most 
probable states given from the density matrix of the sys- 
tem. 

We use the finite-size DMRG algorithm. This consist 
of the following steps: After growing our system up to a 
fixed size L, by means of the infinite-size DMRG, the ba- 
sis of this final system is optimized to best represent the 
desired target state, like the ground state, by sweeping 
through the system repeatedly. A sweep over the system 
is an iterative process which starts with a small block on 
the right extreme of the chain. This is grown to in the 
left direction by adding a site to the right block and con- 
necting it to a bath or environment on the left side. The 
environment information was collected from the infinite- 
size algorithm. The total size of the system is always 
kept constant. As soon as the decreasing size of the left 
block reaches a single site the procedure is stopped. We 
save the information of the right blocks and can use it 
now to start a similar procedure with a block on the left 
side of the chain being grown in the right direction. This 
procedure is repeated until convergence is reached. 

With each step, the chain grows one site in the current 



direction, and the basis of the new system must be trun- 
cated to keep the Hilbert space manageable. All the nec- 
essary operators are transformed and stored every time 
this happens. With every step, the choice of states in the 
truncation of the basis becomes a better representation 
of the system. This leads to an optimal truncated basis 
for representing the target state on the finite system. Af- 
ter convergence was reached, we can proceed to measure 
other observables. 

The numerical error caused by truncation of the orig- 
inal basis can be measured through the weight of the 
states that were discarded in a DMRG step. Our sys- 
tems, with L = 240 sites under open boundary condi- 
tions, were investigated keeping m = 256 density-matrix 
states, rendering a maximum truncation error of approx. 
1(T 6 . 



Density-density correlation function 

An operator A, acting either on the left or on the 
right block, can be written in the basis of the specific 
block as (^l-AI 1 ^). In the case of correlation functions 
like ($|AB|^), handling operators requires some extra 
attention. The operators A and B can operate either on 
equal or on different blocks. The last case may lead to 
errors in the calculation of the expectation value of the 
product AB, since each operator is separately written in 
its corresponding basis. The way to proceed is to build 
the exact operator C = AB, in a full basis from the be- 
ginning and transform it as is done for the rest of the 
operators. 

We calculated the TL parameter K p by measuring 
the correlation function between the sites x and Xq: 
C x = (n(x)n(xo)} — (n(x)){n(xo)}, where the static ex- 
pectation values were subtracted. To reduce the effect 
of the local density oscillations, we take the average over 
pairs of correlation functions for neighboring sites calcu- 
lating C(r) = (C x + C x +i)/2, with r = \x — Xo\ and xq 
in the middle point of the chain. Due to the symmetry 
of the problem we can, in principle, choose either branch 
of the system to estimate K p . 



IV. RESULTS 

Using systems with open boundary conditions, finite 
size effects are induced. One example of these effects are 
the local density oscillations and the charge accumulation 
close to the edges of the system, shown in Fig. [3] The 
charge distribution is expected to be symmetric around 
the middle of the chain. We observe, however, that the 
symmetry is slightly perturbed, as seen in Fig. [31 at the 
positions where the Coulomb potential switches values. 
For our purposes, such small changes are negligible, spe- 
cially after taking the average over pairs of correlation 
functions, as explained in Sec. IIII1 It is still observed 



that the charge density remains fairly homogeneous in 
the valley of the Coulomb interaction. 
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FIG. 3: Density profile (n{x)) for Hetero structure I, where 
the on-site Coulomb potential (bottom line on the graph with 
scale on the right) is Ul — Ur — 1.1 and Uc = 0.9. V(x) — 
0.0 for all sites. The band filling is n = 0.5. 

To estimate K p , we fit the values of the numerical data 
to the leading term of Eq. @, K p /TT 2 r 2 , leaving aside, 
at least for the moment, the logarithmic corrections. In 
figures[5]and[6]the density-density correlations are shown 
for both heterostructures and for band fillings of n = 
0.20 and n — 0.50. We will refer to the region from the 
middle of the chain up to the boundary (at the site xr) 
of the Coulomb valley as the region R%, and from this 
point until the end of the chain as the region i? 2 - I n 
the following we describe with detail the results for each 
heterostructure. 

Density-density correlation functions for Het- 
erostructure I. As a first structure we take a slightly 
inhomogeneous Hubbard lattice setting Ul = Ur = 1.1 
and Uc = 0.9 with V{x) = 0.0 for all sites. The valley 
in the on-site Coulomb repulsion has sharp edges at the 
sites Xl and xr, as shown in the Fig. [3J This however, 
and as seen from the full line in both Fig. [5] and in Fig. [HI 
does not alter significantly the continuous decay of the 
correlation function. For all band fillings, 0.10 < n < 1.0, 
the power-law decay extends beyond the boundary point 
and is not completely constrained to any of the regions 
Ri or i?2- In Fig. [7] the values for the TL parameter 
are shown as a function of the band filling. We observe 
that K p < 1.0, which indicates that spin or charge den- 
sity waves are present. The 2fcp oscillations can be also 
observed in the graphs and a closer view is presented in 
Fig. 2J A fitting of the 2fcp oscillations succeeded over the 
whole system only for n < 0.5. In the case of n > 0.5, 
the fitting of the data was only successful at large dis- 
tances. This behavior is reflected on the values of K p , as 
we observe in Fig. [7] the two different values sets for the 
density intervals already mentioned. With this we con- 
firmed the power-law decay of the correlation functions 



as it was possible to determine K p also including the first 
logarithmic correction. 

We compared the results with a similar configuration, 
this time with an interaction of the form U{x) = cos(ax) 
with a a constant. The valley around the center of the 
system remains but the transition on the potential to- 
wards the ends is done in a smoother way. This varia- 
tion of U resulted in the same values for the correlation 
functions as already presented, showing that the sharp 
edges of the on-site potential did not influence strongly 
the Luttinger liquid behavior of the system. 
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FIG. 4: Density-density correlation function for Heterostruc- 
ture I with n — 0.50. The crosses show the numerical 
data and the solid line result after fitting including the term 

A 1 cos(2k F x)x- {1+K " ) ln(x)~i. 
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FIG. 5: Density-density correlation function for Heterostruc- 
tures I (continuous line) and 77 (broken line). In the first case 
we found K p — 0.864. In the second case, K p = 1.093 only 
within Ri, with n — 0.20 in both cases. 

Density-density correlation functions for Het- 
erostructure II. In this case we keep the valley in the 
Coulomb interaction of the former case: Ul — Ur = 1.1 
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FIG. 6: Density-density correlation function for Heterostruc- 
tures I (continuous line) and 77 (broken line). In the first case 
we found K p = 0.838. In the second case K p — 1.15 only 
within Ri, with n — 0.50 in both cases. 



and Uc = 0.9. Furthermore we simulate two poten- 
tial walls by introducing the confining potential V XL — 
V XR = 10.0 and V(x) = 0.0 for the rest of the sites. 
The results in the case of the Heterostructure II distin- 
guished strongly from those previously described. The 
introduction of the confining potential V(xl) and V(xr) 
generated stronger changes from one region to the other, 
killing the oscillations beyond the xr point. In Fig.[5]and 
Fig. [5] the broken line shows how the decay of the corre- 
lation function is abruptly interrupted by the scattering 
potential, not having further space to fully establish the 
decay in the amplitude of the 2A:p oscillations. We found 
however, that even in the constrained region R\ the cor- 
relation function obeys a power law decay with K p > 1.0, 
which is indicative of a Fermi liquid behavior. 



V. CONCLUSIONS 

In this paper we have investigated the behavior of den- 
sity correlation functions in one-dimensional heterostruc- 
tures. We described how junctions between different 
types of atoms influence the variation in space of the 
TL parameter. The heterostructures as defined can be 
seen as unions of subunits with different coupling con- 
stants in which the TLL for homogeneous systems is to 



be expected. However, our findings show that a slow 
variation of the on-site Coulomb potential, as in the first 
case, does not interrupt nor split the decay of the density- 
density correlation functions between the regions and the 
system as a whole behaves as a TLL. Similar systems 
were investigated [16| where the on-site Coulomb poten- 
tial was turned on and off over the subchains. For such 
systems an effective exponent, K* = f(Ki >p ,K2 tP ), was 
calculated considering, for example, two subchains which 
were assumed as independent, homogeneous TLL's. Us- 
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FIG. 7: TL parameter K p for both heterostructures as a func- 
tion of the band filling. For Heterostructure II , only within 
the region R\. For Heterostructure I K p < 1.0 indicating 
a Luttinger liquid behavior. In the second heterostructure 
K p > 1.0, which is the benchmark of a Fermi liquid. 



ing DMRG, their reported values could only be partially 
reproduced, namely for densities n < 0.6. 

We have also found a completely different behavior re- 
sulting from the introduction of a scattering potential V 
at the junctions between the subunits, as done for the 
systems in the second case. Our findings in such case 
show that the TLL is not a universal feature for one- 
dimensional systems. Concerning the dynamics in het- 
erostructures, further work remains to be done. Trans- 
port properties at temperatures different from zero are 
a key in the construction of properly tunable electronic 
devices. 

We wish to acknowledge useful discussions with A. Mil- 
lis and thank C. Kollath for a helpful revision of the pa- 
per. 
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